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Preface 


This  thesis  represents  the  result  of  a ten  month  inves- 
tigation of  a second  order  technique  for  finding  the  optimum 
controls  to  be  used  by  an  aircraft  evading  an  air-to-air  heat 
seeking  missile.  The  approach  involves  the  application  of 
optimal  control  theory  utilizing  a Differential  Dynamic  Pro- 
gramming Model  to  determine  the  optimum  controls. 

This  thesis  was  sponsored  by  the  Aeronautical  Systems 
Division  as  part  of  a study  being  conducted  by  the  Air  Force 
for  a better  visual  display  to  the  pilot  for  evading  an  air- 
to-air  missile . 

I have  had  great  satisfaction  in  developing  the  algo- 
rithm for  this  project.  My  only  disappointment  is  that  I 
was  unable  to  finish  the  algorithm  in  order  to  obtain  the 
desired  results;  however,  I do  feel  that  I have  laid  the 
groundwork  that  could  lead  to  a worthwhile  second  order 
technique . 

I wish  to  sincerely  thank  my  thesis  advisor,  Major  James 
Funk,  for  his  assistaiice  and  guidance  during  this  project.  I 
would  also  like  to  thank  my  sponsor,  Mr.  Mike  Breza,  for  his 
helpful  suggestions. 

I dedicate  this  thesis  to  my  wife,  Sharon,  and  my  chil- 
dren, Todd  and  Kristi,  who  gave  me  encouragement  and  under- 
standing, and  exercised  an  incredible  amount  of  patience. 

Robert  Smith 
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Abstract 


The  purpose  of  the  study  is  to  formulate  a method  to 
determine  the  control  strategies  that  maximize  the  probabil- 
ity of  sui^vival  for  an  evading  aircraft.  This  is  equivalent 
to  minimizing  the  probability  of  kill  for  the  attacking  air- 
to-air  missile.  The  controls  of  the  evading  aircraft  consist 
of  the  commanded  angle  of  attack,  bank  angle,  and  the  com- 
manded coefficient  of  tirust.  The  missile  model  developed 
is  a typical  air-to-air  infrared  missile  using  proportional 
navigation  steering.  The  probability  of  kill  is  modeled  as 
ellipsoidal  iso-cost  surfaces  with  a cost  value  that  decays 
exponentially  as  the  ellispoid  size  increases.  The  flatten- 
ing and  position  of  the  ellipsoid  centroid  account  for  the 
shape,  orientation  and  ‘/ulnerability  of  the  aircraft.  The 
problem  terminates  when  the  line-of-sight  from  the  missile  to 
the  target  aligns  with  the  surface  of  the  missile’s  fuzing 
cone.  The  algorithm  developed  employs  a second  order  dif- 
ferential dynamic  programming  model  for  optimizing  the  con- 
trols of  the  evading  aircraft. 
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OPTIMAL  MISSILE  EVASION 


I.  Introduction 


Purpose  of  the  Study 

With  the  newer  fighter  aircraft  and  fire  control  system 
a problem  exists  for  both  a+  eking  and  evading  aircraft.  A 
method  of  determining  and  iisplaying  real  time  probability  of 
kill  (Pj^.)  to  the  pilot  of  the  attacking  aircraft  using  an  air- 
to-air  missile  is  very  desirable.  On  the  other  hand,  the  ef- 
fectiveness of  missiles  against  aircraft  has  made  it  extremely 
desirable  to  provide  the  pursued  pilot  with  aids  for  evading 
missiles.  This  study  will  be  concerned  with  the  latter  of 
these  two  problems,  that  is  the  problem  associated  wdth  the 
e-'  ling  aircraft. 

The  evasion  problem  is  a two  step  process  consisting  of: 

1,  Determining  the  optimum  control  strategies  to 
be  used  by  the  evading  pilot. 

2.  Providing  on-board  computation  and  display  of 
cues  or  solutions  for  the  pilot. 

This  study  represents  the  first  step  in  that  process,  that  is 
determining  the  optimum  controls  to  be  used  by  the  evading 
pilot. 

Background 

The  maximum  and  minimum  effective  ranges  of  air-to-air 
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missiles  are  functions  of  attacker  state,  target  state,  and 
performance  characteristics  of  both  vehicles.  Operational 
experience  has  demonstrated  that  pilots  have  difficulty  in 
accurately  estimating  valid  launch  conditions  during  an  en- 
counter. Consequently,  effective  employment  of  air-to-air 
missiles  requires  two  distinct  functions;  1)  an  accurate 
missile  launch  envelope  computation  performed  in  an  airborne 
computer,  and  2)  display  of  appropriate  parameters  to  the 
attacking  pilot  so  that  he  can  recognize  and  take  advantage 
of  valid  launch  opportunities. 

The  missile-launch  envelope  considerations  are  really 
beyond  the  scope  of  this  study,  which  concentrates  on  eva- 
sion, However,  the  solutions  to  the  evasion  problem  could 
be  used  as  a basis  to  determine  launch  envelopes.  It  is 
still  a sizeable  step  to  solve  and  process  the  necessary 
solution  data  in  order  to  obtain  usable  prcbabj.lity  of  kill 
information. 

Present  airborne  digital  computers,  such  as  in  the  F-15, 
have  made  it  practical  to  develop  and  to  implement  evading 
strategy  computations  and  to  display  these  parameters  to  the 
pilot.  Thus  this  study  is  concerned  with  the  development  of 
improved  strategies  to  be  used  by  an  evading  target;  solu- 
tions that  consider  more  informat'.on  than  past  methods. 

Scope 

The  state  variable  equations  are  simplified  where  ever 
possible  without  sacrificing  any  significant  realism.  There 
is  not  as  much  freedom  of  motion  as  in  the  actual  case  due  to 


the  assumption  of  coordinated  turns  by  the  aircraft. 

The  set  of  controls  were  chosen  such  that  the  evading 
aircraft  has  flexibility  in  controls  and  is  still  realistic. 

The  terminal  cost  function  corresponds  generally  to  a 
detailed  simulation  of  the  end  game  which  incorporates  the 
shape,  orientation  and  vulnerability  of  the  aircraft.  Indi- 
vidual component  vulnerabilities  of  a particular  aircraft  are 
not  considered,  only  the  general  characteristics  of  a typical 
aircraft. 

The  second  order  differential  dynamic  programming  algo- 
rithm was  selected  for  evaluation  in  solving  this  problem. 
This  choice  was  based  on  the  high  degrtc  of  nonlinearity  of 
the  state  equations,  and  convergence  difficulties  with  a 
first  order  algorithm.  The  algorithm  was  adapted  to  this 
problem  using  unspecified  final  time. 

This  report  will  be  limited  to  an  approach  for  obtain- 
ing the  following  informations 

1.  The  trajectory  of  the  missile. 

2.  The  optimum  trajectory  of  the  evading  target. 

3.  The  minimum  Pj^  of  the  missile. 

4.  The  optimum  controls  used  by  the  evading  air- 
craft . 

The  information  obtained  is  determined  from  a specific  set  of 
launch  conditions. 

Assiimptions 

Certain  assumptions  can  be  made  to  reduce  *he  complexity 
of  the  problem  without  significantly  affecting  the  character 
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of  tht  solution.  Thoy  are  as  follows t 


1.  Rigid  body  dynsmio  moaels  for  both  die  missile 
and  the  evading  target. 

^2.  The  yaw  angle  is  assumed  negligible. 

3.  For  the  evader,  the  angle  of  attack  and  thrust 
respor.i’^3  to  respective  commands  are  modeled 
as  lineal  first  order  systems  with  api:iropriate 
time  oov.stants. 

4.  The  bank>i.'rv-,i6  tine  responses  are  assumed  rap- 
id eno'«:.9,h  to  neglect  any  time  delay  in  the  re- 
sponse . 

5.  The  missile  fusing  cone  angle  is  fixed,  as  is 

the  delay  time. 

C,  The  densi'i  jT  tne  atmosphere  as  a function  of 
altitude  is  &iven  byi 


where 

pQ  * .0023769  slugs/ft^ 

* 23800  ft 

Z " the  altitude  above  the  earth's  sur- 
face. 


General  Approach 

The  material  is  presented  in  the  following  order.  First 
the  equations  of  motion  for  both  the  missile  and  evading  air- 
craft are  derived.  The  proportional  navigation  steering  is 
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then  developed  along  with  the  computed  acceleration  for  the 
missile.  The  state  equations  for  the  dynamic  model  are  next 
outlined.  This  is  followed  by  the  derivation  of  the  cost 
function  derived  from  the  probability  of  kill  geometry.  Next 
the  terminal  constraint  geometry  is  obtained.  The  required 
equations  for  the  differential  dynamic  programming  algorithm 
(Ref  1i47)  are  outlined.  Finally,  the  computational  proce- 
dure used  for  the  algorithm  is  discussed,  followed  by  results 
and  conclusions. 
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The  State  Eciuations 


General  Description 

The  continuous -time  dynamic  system  modeled  for  this  pro- 
gram is  described  by  the  following  set  of  nonlinear  ordinary 
differential  equations  <' 

= f(X^,u)  I X^(tQ)  = X^^ 

\ ” \ 

where  the  subscripts  T and  m denote  target  states  and  missile 
states  respectively.  The  perfonnance  of  the  system  is  meas- 
ured by  minimization  of  a terminal  cost  function  given  ast 

s:  (2) 

subject  to  the  terminal  constraint  of  the  formi 

i|;(X^(tj),X^(t^))  = 0 (3) 

where  the  final  time  t^  is  given  implicitly, 

Defining  the  State  Vector 

The  state  vector  for  this  optimal  missile  evasion  prob- 
lem includes  the  distance  components  between  the  target  and 
missile,  the  velocity  of  both  target  and  missile,  the  heading 
of  target  and  missile,  the  flight  path  angle  of  target  and 
missile,  the  angle  of  attack  of  the  target,  and  the  coeffi- 
cient of  thrust  of  the  target.  In  standsird  notation  the 
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bodies  with  their  velocity  vectors  expressed  in  an  inertial 
reference  frame  (ijk)  as  depicted  in  Pig.  1. 


Fig,  1.  Velocity  Vector  in  Inertial  Frame 

The  velocity  axis  frame  is  designated  as  e^,  e^.  The 
velocity  of  the  center  of  mass  in  the  velocity  axis  system  is 

V = (6) 

transformation  into  the  inertial  reference  frame  yields 

X = Vcos(y)  sin(i|;) 

Y = Vcos(y)  cos(i|;)  (7) 

Z = Vsin(Y) 

where  X and  Y are  horizontal  position  coordinates,  Z is  alti- 
tude, V is  speed,  y is  the  flight  path  angle  and  ij;  is  the 


heading  angle  with  respect  to  the  Y axis. 

In  order  to  obtain  the  accelerations  of  each  vehicle, 
the  derivative  of  the  velocity  in  the  ,floving  frame  must  be 
taken. 

^ Ve^  (8) 

and  recalling  that 

= (3  X e^  (9) 

v^ere  the  angular  rate  co  is  given  as 

55  « Ye,|,  - (10) 


and  ^ere 


e^  a cos(Y)e^  + sin(Y)eY  (11) 

therefore,  combining  Eqs  (9).  (10),  and  (11)  and  substituting 
into  Eq  (8)  yields 

^ = Ve^  + Vye^  + ViircoB(Y)e^  (12) 

The  acceleration  of  the  body  times  the  mass  of  the  cody  is 
equal  to  the  sum  of  the  external  forces,  or 

= f (13) 

where  f is  composed  of  thrust,  lift,  drag,  and  gravity.  A 
more  detailed  representation  of  the  velocity  axis  coordinate 
system  with  the  forces  acti.ig  on  the  vehicle  is  shown  in 
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Fig.  2 amd  Fig.  3* 


aaangle  of 

attack 

Da drag 

Chord  Line 

g=gravity 

Lalift 

T= thrust 

Fig.  2.  Velocity  Axis  Coordinate  System 
With  Angle  of  Attack 


3=bank  angle 


Pig.  3.  Velocity  Axis  Coordinate  System 
With  Bank  Angle 


Referring  to  Pig.  2 and  Pig.  3 the  following  force  components 
are  obtained i 

Drag  = -D  e^ 

^ = T cos(a)ey  + T sinCa)?^ 

Lift  = L co8(0)e^  + L sin(0)e^ 

Gravity  = -g  cosCy)©^  - S sin(Y)e^  (1^) 

where 

D = 1/2  p S [Cd^  + K(CLjj^  a)^] 

L = 1/2  p S CL^  a (15) 

Combining  Eqs  (12)  through  (15)  the  following  acceleration 
terms  are  obtained  in  the  inertial  axis  coordinate  system t 

Y « (T  C08(a)  - 1/2  p S V^(Cdp  + K(CLj^  a)^))/m 

- g sin(Y) 

Y = 1/2  p S V(CLq^  a)  cosO)/m  + 

- g COS(y)A 

4r  » 1/2  p S Y(CL^  a)  sin(3)/m  cos(y)  (16) 

Required  Aerodynamic  Acceleration  of  the  Missile 

Proportional  navigation  provides  a rate  of  change  of 
the  missile  heading  directly  proportional  to  the  rate  of  rota^ 
tion  of  the  line -of- sight  from  the  missile  to  the  target 
(providing  the  rates  are, within  missile  performance  limita- 
tions). The  line-of-sight  rate  from  the  missile  to  the 
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spum"! 


target  is  defined  as 


je_  = — TT-^ 


(17) 


2 

where  r = r • r and  r is  the  relative  distance  between  the 
missile  and  target,  therefore » 

r = (X^  - X„)I  + (Yj  - X„)l  + (Zj  - Z„)lt  (18) 


and  Vj,  is  the  relative  velocity  between  the  missile  and  tar- 
get, therefore  1 

= (V^  cosCy^,)  sin(i};^)  - ooB(y^)  sinCijij^))! 

+ (Vj  cos(Y|p)  cos(i{;j)  - Vj^  cos(Yjn)  cos((lfj^))1 

+ (V^  sin(Yj.)  - \ sin(Yjn))ic  (19) 


The  desired  turn  rate  for  the  missile  is  then  equal  to  a pro- 
portional navigation  constant  (n)  times  the  line-of-sight 
rate,  or  the  desired  turn  rate  is  n For  a typical  air- 
to-air  missile  n is  in  the  range  of  2 to  and  for  this 
problem  n is  set  equal  to  3.  The  computed  acceleration  of 
the  missile,  defined  as  a^,  is  then  given  asi 


5c  = 


^ Ir 


X V 


m 


(20) 


In  order  to  determine  the  required  aerodynamic  acceleration 
of  the  missile,  designated  as  a^^,  the  effect  of  gravity  must 
be  incorporated.  Therefore,  combining  Eq  (20)  with  the  ef- 
fect of  gravity  the  required  aerodynamic  acceleration  becomes 


.laiiimaiiiiujiMM 


where  is  the  required  aerodynamic  acceleration  and  is 
the  component  of  gravity  normal  to  the  velocity  of  the  mis- 
sile. a^^  is  the  aerodynamic  acceleration  which  the  vehicle 
should  produce  in  order,  to  obtain  the  proper  normal  accelera- 
tion in  the  presence  of  gravity.  Referring  to  Fig.  4,  grav- 
ity in  the  velocity  frame  system  is  given  ast 

g velocity  = -g  sin(Y)e„  - g cos(y)  + Oe.  (22) 
frame  ^ Y 


In  the  required  aerodynamic  acceleration  formula,  the 
effective  component  of  gravity  is  that  component  which  is 
normal  to  the  velocity  vector i or  from  Eq  (22)  effective 
component  of  gravity  becomes i 

g velocity  = -g  cos(Y)e  (23) 

frame  ^ 

effective 

Referring  to  Fig.  1,  transformation  from  the  velocity 
coordinate  system  to  the  inertial  coordinate  system  is  given 
by  I 


» m 

1 

o 

o 

rl. 

m M 

cos(i};)  -sin(ijj)  0 

i 

= 

0 cos(y)  sin(Y) 

sin(i|;)  cob(iIj)  0 

j 

m m 

0 -sln(Y)  cos(y) 

0 0 1 

mm  ■■ 

k 

m m 

which  simplifies  toi 


m m 

@1 

'If 

cos(i}f)  -sin(ijj)  0 

m m 

9 

1 

®v 

s 

cos(y)  sin(i|;)  cos(y)  cos(ij;)  sin(Y) 

» 

J 

e 

Y 

Ml  m 

-sin(Y)  sin(il;)  -sin(Y)  cos(ijj)  cos(y) 

m»  m 

k 

m m 

Combining  Eqs  (23)  and  (24)  the  effective  component  of  gravi- 
ty in  the  inertial  coordinate  system  then  becomes j 

e inertial  = = g cos(y)  sin(Y)  3in(ili)i 

frame  ” ^ 

effective  + g cos(y)  sin(Y)  cos(il;)j 

- g cos^(Y)h  (25) 

Therefore,  Eq  (25)  along  with  the  computed  acceleration  of 
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the  missile  (£q  20)  determines  the  required  aerodynamic  accel- 
eration of  the  missile. 

Angle  of  Attack  of  the  Missile 

A detailed  representation  of  the  forces  acting  on  the 
missile  Is  depicted  In  Fig.  2.  Equating  forces  and  noting 
that  the  thrust  vector  T acts  along  the  vehicle  axis  of  sym- 
metry, the  force  normal  to  the  flight  path  Is  given  ast 

^normal  * ® ^ 

Assuming  that  the  angle  of  attack  remains  small  so  that  sln(a) 
approximately  equals  a,  and  substituting  In  for  lift,  £q  (15), 
the  force  normal  to  the  flight  path  becomes t 

m * 1/2  p S a + T a (27) 

Solving  £q  (27)  for  the  angle  of  attack  of  the  missile 
yields 

a^  m 

« * ^ (28) 

^ 1/2  p V2  S CL,  a + T 

d 

Bank  Angle  of  the  Missile 

In  order  to  determine  the  bank  angle  of  the  missile, 

refer  to  Fig.  5 for  a schematic  of  the  bank  angle  In  the 

velocity  coordinate  system.  The  missile  bank  angle,  3^,  Is 

the  angle  between  the  unit  vector  e and  the  vector  u«  where 

Y ^ 

is  the  unit  vector  along  the  direction  of  the  required 
aerodynamic  acceleration  a^.  Using  the  law  of  cosines,  the 
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Pig,  5*  Determination  of  Missile  Bank  Angle 


bank  angle  is  then  determined  byj 


= ± cos“^(e^  • \) 


>diere  the  unit  vector  is  determined  byj 

3.  = A- 

“ l5„l 


(29) 


(30) 


The  -ttiird  component  of  Eq  (24)  gives  e^ 

e^  s -Bin(Y)  sin(i|r)i  - sin(Y)  co8(i|;)q  + cos(Y)ic  (31) 

Referring  to  Pig,  5 and  Eq  (29),  the  following  stipulations 
are  made  on  the  bank  angle t 

1,  If  e^  • Ug^  > 0 then  is  between  0 and  180 
degrees, 


2.  ® then  3j^  is  between  180  and  360 

degrees, 


v^ere  the  unit  vector  e,^  is  given  ast 

e,  = cos(t|;)I  - sin(!j;)^ 


(32) 


In  other  words,  the  positive  sign  of  is  used  when  condi- 
tion 1.  is  satisfied,  and  the  negative  sign  of  3^  is  used 
when  condition  2.  is  satisfied. 


the  Nonlinear  Differential  State  Equations 


The  velocity  equations  represented  by  Eq  (7)  and  the 
acceleration  equations  represented  by  Eq  (16)  produce  the 
following  set  of  nonlinear  ordinary  differential  state  equa 
tions  used  in  this  program t 

e • 

Xj,  - Xj  = Vjg  cos(Yja)  8in(i|;gji  - cos(yt)  sin(!jfj) 


\ ’ ^T 


= Vj^  cos(Ynj)  ■ ^T  cos(i|fj) 


~ 8in(Yij«) 


= (T„  oos(a„)  - 1/2  p„  v/  S„(0a„„ 
♦ - 6 
- 1/2  pj  Sj(0j  co8(aj)  - (Od^j 

+ KjCOLjjj  aj)*))/mj  - g 3in(Yj) 

= Pm  ®m  \ Ol'am  “n  »i"<8m>/'V 


Ym 


“ 1/2  Pu  CL^ip  Qtjj,  8ln(3ip)/ni^  cos(y^) 

“ V2  P„  S„  V„ 

* ’'n.  % - 8 ‘=°“<lf|n>/''n. 

■“  P«p  S(p  cos ( 8(p )(i(p/in,j) 

+ C^  sin(a^)/m^)  - g cos{y^)/y^ 


a. 


T 


(ot^  “ Clip)/Tp 

c 


(Cq,  •*  C(p)/T(ji 
C 


(33) 


The  subscripts  T and  m denote  the  target  and  missile  parame- 
ters respectively,  Tp  is  the  time  constant  for  the  angle  of 
attack  response,  and  is  the  time  constant  for  the  thrust 
response  of  the  target,  ajj^  and  have  previously  been  de- 
fined by  Eqs  (26)  and  (29),  , 0^,  and  Cip  are  the  control 

O 0 

variables  defined  by  Eq  (5).  and  where  and  p^  are  defined 
by  I 


Pm  “ Po® 


= p e(-V^'o) 


(34) 
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III.  Terminal  Cost 


The  objective  of  the  terminal  cost  function  is  to  make 
the  model  correspond  to  detailed  simulation  of  the  end  game 
which  incorporates  the  shape,  orientation  and  vulnerability 
of  the  aircraft.  Individual  vulnerability  of  a particular 
aircraft  is  not  considered,  but  only  the  generic  character. 
The  ellipsoid  parameters  were  chosen  to  fit  detailed  sim- 
ulation results  reasonably  well,  The  general  form  of  the 
terminal  cost  is  given  byj 

(2) 

As  previously  stated,  this  optimisation  problem  requires 
the  minimization  oi  this  terminal  cost.  The  terminal  cost  is 
a convex  function  described  by  ellipsoidal  constant  cost  sur- 
faces centered  about  a reference  point  near  or  on  the  air- 
craft. Each  surface  represents  a constant  value  for  the 
probability  of  kill  for  a particular  missile.  The  probabil- 
ity of  kill  (Pj^)  values  decrease  exponentially  with  increas- 
ing concentric  ellipsoid  size.  The  ellipsoids  can  be  con- 
sidered fixed  with  respect  to  the  target  vehicle  - moving  and 
rotating  with  it.  The  Pj^  for  each  missile  approach  path  is 
a function  of  the  crossing  aspect  angles  and  the  separation 
distance  between  the  missile  and  target  paths,  as  well  as  the 
relative  velocity  and  altitude. 

The  Pj^  for  each  flight  is  evaluated  by  determining  the 
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separation  distance  of  the  missile  flight  line  projection 
from  the  target  at  the  terminal  time.  The  separation  dis- 
tance of  the  missile  flight  line  projection  is  designated  r 
in  Pig.  6^ and  will  be  more  precisely  defined  in  the  following 
paragraphs. 


Pig.  6.  Terminal  Cost  Geometry 


Designating  the  fuzing  time  of  the  missile  as  t^,  which 
is  the  time  when  the  target  intercepts  the  missile  fuzing 
cone  and  is  the  beginning  of  the  missile  delay  interval,  T^, 
defined  as  the  period  of  time  from  fuzing  to  detonation.  The 
target  position  at  detonation,  designated  as  the  origin  0, 
can  be  determined  by 

0 = X,j,(t^)  + V,j,(tjp)  • T^  (35) 

where  X,p  is  the  target  position  at  t^  and  V,p  is  the  target 
velocity  at  t^,  in  the  inertial  coordinate  system.  The  sep- 
aration distance  r is  the  offset  distance  of  the  missile 
flight  path  projection  from  the  origin.  The  vector  P,  which 
is  the  line-of-sight  vector  from  the  missile  at  fuzing  to  the 
target  at  detonation  is  then  given  by: 

P = - 0 (36) 

where  X^^^Ct^)  is  the  position  of  the  missile  at  the  fuzing 
time.  Substituting  Eq  (35)  into  Eq  (36)  yields 

P = iTi(tf)  - - V,j,(t^)  • Tf  (37) 

The  vector  r,  which  is  the  vector  from  the  origin  per- 
pendicular to  the  missile  flight  line  projection,  is  a func- 
tion of  P, 

r = (u  X P)  X u (38) 

which  is  equivalent  tos 
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Designating  the  fuzing  time  of  the  missile  as  t^.,  which 
is  the  time  when  the  target  intercepts  the  missile  fuzing 
cone  and  is  the  beginning  of  the  missile  delay  interval, 
defined  as  the  period  of  time  from  fuzing  to  detonation,  The 
target  position  at  detonation,  designated  as  the  origin  0, 
can  be  determined  by 

0 = X^(tj.)  + V^(t^)  • T^  (35) 

where  is  the  target  position  at  t^  and  V|p  is  the  target 
velocity  at  t^,  in  the  inertial  coordinate  system.  The  sep- 
aration distance  r is  the  offset  distance  of  the  missile 
flight  path  projection  from  the  origin.  The  vector  P,  which 
is  the  line-of-sight  vector  from  the  missile  at  fuzing  to  the 
target  at  detonation  is  then  given  byt 

P = X^(t^)  - 0 (36) 

where  X^('tj)  is  the  position  of  the  missile  at  the  fuzing 
time.  Substituting  Eq,  (35)  into  Eq  (36)  yields 

I = ^m^^f^  - ^T^^f^  “ ^T^^f^  * "^f 

The  vector  r,  which  is  the  vector  from  the  origin  per- 
pendicular to  the  missile  flight  line  projection,  is  a func- 
tion of  P, 

r = (u  X P)  X u (38) 

which  is  equivalent  toi 
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(39) 


r = P - (P  • u)  u 

nAiere  u is 
following i 


Combining  Eqs  (37) i (39).  and  (40)  yieldsi 

E ' « - Tf  Vt  - (iX  • }!«)  \/\^ 

* TfA„^  C^i) 

where  all  values  are  determined  at  the  fuzing  time  t^,  and 
where  • designates  the  vector  dot  product. 

The  ellipsoidal  equation  for  iso  - Pj^  contours  in  matrix 
form  ist 

* r®  F r (42) 

where  r,  in  Eq  (4l),  can  be  expressed  in  the  inertial  coordi- 
nate system.  F is  most  easily  defined  in  the  target  coordi- 


nate  system 

0 0 

Ft  - 

0 l/dy2  0 

0 0 

where  F is  the  scaling  matrix  determined  to  approximate  the 
Pjj.  data  of  detailed  vulnerability  studies  by  the  ellipsoid 


the  missile  unit  velocity  vector  determined  by  the 


u « 


(40) 
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surfaces.  In  order  to  express  F in  the  inertial  coordinate 
system,  the  following  similarity  transformation  must  be  made: 


P = oj  Fj  oj 


(W) 


is  the  transformation  matrix  from  the  inertial  coordinate 
system  to  the  target  coordinate  system?  and  transforms 
from  target  to  inertial  coordinates.  In  determining  the 
order  of  rotation  is  through  the  heading  angle  fol- 

lowed by  flight  path  angle  (y^),  and  then  the  bank  angle 
(0^),  and  finally  through  the  angle  of  attack  (a,p) : 


- 

Ct  - 


M * 

O 

O 

I 

cos(0^)  0 -sin(0j) 

0 cos(aj)  sin(aj) 

• 

0 10 

0 -sinCa^p)  cos(a,p) 

sin(3^)  0 cos(3,p) 

o 

o 

I 

cos((|;j)  -sin(il;j)  0 

0 cos(y^)  sinCvip) 

t 

sinCij;^)  cos(i}f^)  0 

0 -sin(Y^)  cosCyj) 

0 0 1 

(45) 


Therefore,  the  transformation  from  the  target  coordinate 
system  to  the  inertial  coordinate  system  is  given  by: 


ol  = (0j)» 


(46) 


Finally,  the  terminal  cost,  defined  as  F(X(t^)?  t^.)  is 
given  by: 

.2 


F(X(t^)?  tf)  = e 


-R‘ 


(4?) 
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2 

where  R is  defined  by  Eq  (42). 

The  type  missile  used  for  this  simulation  is  an  air-to- 
air  heat  seeking  missile.  The  tail  of  the  aircraft  was  chosen 
as  the  ai^  point  (or  origin)  for  the  t2a‘get.  Since  the  mis- 
sile is  guided  towards  the  tail  of  the  evading  aircraft,  kill 
probabilities  are  higher  if  the  detonation  is  forward  rather 
than  behind  the  target  aim  point.  A bias  term  is  used  to 
shift  the  center  of  the  ellipsoids  forward  to  account  for 
more  probable  fuzing  forward  of  the  aim  point.  The  amount 
of  shift  is  designated  as  b,  then  the  vector  from  the  ellip- 
soid center  to  the  tip  of  the  r vector  is 

D = r - b (48) 

where  r is  the  "closest  approach  point"  of  the  projected  mis- 
sile flight  line  projection  to  the  aim  point  (the  analysis 
was  based  on  r as  a parameter).  The  missile  used  for  this 
algorithm  has  a shift  vector  b ofj 

by  = 6 ft 

bj  = 2 ft  (49) 

and,  therefore,  the  ellipsoid  equation  for  the  iso-Pj^  con- 
tours in  matrix  form  become si 

F D (50) 

I 

For  the  missile  used,  the  weighting  matrix,  F^,  in  the 
target  coordinate  system  isi 
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F = 


1 

(21)' 


(22) 


(16)' 


(51) 


Therefore,  the  terminal  cost  for  the  missile  used  in  this 
program  hecomest 


(2i) 


-D^CJ 


1 

(22)‘ 


1 

(16)‘ 


F(X(t^)l  t^)  * 


C|  D 


(52) 


4 ^ 


IV,  Terminal  Constraint 


A typical  air-to-air  missile  has  a fuzing  cone  angle 
(PCA)  of  approximately  6o  degrees.  The  proportional  naviga- 
tion steering  attempts  to  maneuver  the  missile  so  that  the 
target  is  constrained  inside  the  fuzing  cone  angle  at  all 
times.  During  the  terminal  portion  of  the  flight,  as  the 
missile  approaches  the  target,  practical  limitations  on  mis- 
sile maneuvering  will  allow  the  target  to  reach  the  fuzing 
cone  angle.  At  the  time  when  the  line-of-sight  from  the  mis- 
sile to  the  target  lies  on  the  missile  fuzing  cone  the  missile 
‘fuzing  delay  is  initiated,  and,  following  a preset  delay  time, 
detonation  is  programmed  to  occur.  The  fuzing  cone  angle  and 
delay  times  are  normally  chosen  to  give  "good”  fragment  pat- 
terns . 

The  terminal  constraint  for  this  problem  is,  therefore, 
when  the  line-of-sight  from  the  missile  to  the  target  equals 
the  fuzing  cone  angle,  Refer  to  Fig,  7 for  a schematic  of  the 
terminal  constraint. 

Let  A equal  the  unit  missile  axis  vector  and  S equal 
the  line-of-sight  (LOS)  vector  from  the  missile  to  the  target} 
therefore,  when  the  LOS  lies  on  the  edge  of  the  fuzing  cone 
angle  then 


S X A = |S|  |A|  sin(PCA)  (53) 

where  lA|  is  equal  to  1.  S which  is  the  LOS  vector  from  the 
missile  to  the  target  is  then  equal  to 
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Pig.  7,  Terminal  Constraint  Geometry 


(5^) 


The  terminal  constraint,  designated  as  t^),  is  then 

given  as I 

i!;(X(tj),  t^)  = Is  X A|^  -ISl^sin^CFCA)  (55) 


A which  is  the  unit  missile  axis  vector  in  the  missile  coordi- 
nate system  is  equal  to 


A = 


0 

1 

0 


(56) 


In  order  to  obtain  A in  the  inertial  coordinate  system  the 
transformation  from  the  missile  coordinates  to  the  inertial 
coordinate  system  must  be  made.  The  transformation  from  the 
missile  to  the  inertial  coordinate  system,  Aj,  is  of  the  same 
form  as  the  transformation  matrix  for  the  target  Aj,  and  re- 
calling that 

A™  = (aJ)*^  (57) 


then  Aj  becomes 
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1 

o 

o 

cos{0j|^)  0 -sin(3j^) 

0 cos(aj^)  sin(aj^) 

« 

0 10 

0 -sin(ajjj)  cos(aj^) 

mt  w 

sin(9j^)  0 cos(ej^) 

H*' 

O 

o 

1 

« 

cos(t{jj^)  -sinCilijjj)  0 

■ 

0 cos(Yj^)  sin(Y^) 

• 

sin(il;j^)  cos(i|tj^)  0 

0 -sin(Yj^)  cos(Yin)^ 

0 0 1 

(58) 


where  the  order  of  rotation  is  the  angle  of  attack,  then  bank 
angle,  followed  by  flight  path  angle,  and  finally  by  a head- 
ing change.  Finally,  A in  the  inertial  coordinate  system  is 


(59) 


The  terminal  constraint  equation  is  then  defined  by  Eq  (55) 
where  S and  A are  defined  by  Eqs  (5^)  and  (59)  respectively, 
and  FCA  depends  on  the  missile  being  used.  For  the  missile 
used  in  this  program  the  fuzing  cone  angle  is  60  degrees. 
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Differential  Dynamic  Programming  Equations 

The  second  order  algorithm  for  fixed  end  point  problems 
with  the  final  time  (t^)  given  implicitly  is  discussed  in  Ref 
(It 38).  The  derivation  of  *he  required  equations  will  not  be 
given  here,  but  the  equations  to  be  used  for  this  particular 
problem  will  be  discussed. 

Initially,  a first  order  algorithm  was  considered  for 
use  in  this  program!  therefore,  the  following  sample  problem 
was  attempted  with  the  first  order  method t 

*1*  *2 

= -Xj^  - Xj^3  ♦ u (60) 

The  rate  of  convergence  on  this  nonlinear  problem  using  the 
first  order  method  was  very  slow.  The  second  order  algorithm 
was  then  attempted  on  the  same  problem  with  the  rate  of  con- 
vergence being  greatly  improved  over  that  of  the  first  order 
method.  The  optimal  missile  evasion  problem  is  quite  nonlin- 
ear and  as  previously  outlined  has  a 12  dimensional  state  vec 
tort  therefore,  the  second  order  algorithm  was  selected  even 
though  of  the  disadvantage  of  taking  a large  number  of  second 
order  derivatives. 

The  second  order  method  consists  of  two  phases t the 
first  phase,  optimization  phase,  uses  backward  integra- 
tion of  the  costate  equations  to  determine  the  sensitiv- 
ities of  the  cost  function  for  proper  adjustment  of  the 
controls.  Throughout  this  phase  the  value  of  the  termi- 
nal constraint  functiur.  is  allowed  to  wander.  The  second 
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phase,  which  is  designated  as  the  restoration  phase,  starts  with 
con'trols  obtained  in  the  optimisation  phase  and  then  attempts 
to  restore  the  constraint. 

In  general,  the  function  to  be  minimized  can  be  given 
in  the  general  form  ofi 


4-  P(X(t^)|  t^) 

+ b iIf(X(t^),  tf) 


(61) 


where  the  final  time  t^  is  given  implicitly.  For  this  par- 
ticular problem,  the  controls  are  indirectly  constrained  by 
the  integral  cost  term 

L(X.  u,  t)  = l/2(RA.(aj  + RC-(Cj  )^)  (62) 

c c 

where  RA,  RB,  and  RC,  the  weights  put  on  the  controls,  are 
selected  as  the  inverse  of  the  maximum  value  expected  for 
that  control.  For  this  problem  the  values  are 

RA  =»  10 
RB  = .1 

RC  ^ 73  (63) 

F(X(tj),  tj.),  the  termimd  cost,  is  defined  by  Eq  (52), 

’|f(i('tf,  t^)»  the  terminal  constraint,  is  defined  by  Eq  (55), 
b is  Ihe  time  invariant  Lagrange  multiplier,  and  the  state 
variable  differential  equations  are  of  the  form  X * f(X,  u,  t). 
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' fWli'StlB'W'''  ’"iWPH'n 


The  necessary  sensitivity  equations  to  be  integrated  in  order 
to  minimize  Eq  (61)  are 


H - H(3r,  u.  t) 


Hx  + Vxx(f  - f(X.  Uf  t)) 


^Xb 

= <^x^ 

fu 

''xb 

^xt^ 

= (^x-^ 

fu  »!>’' 

''xt; 

^btf 

= -V 
''xb 

f H 
u uu 

^ f ' 
u 

^bb 

= -V 
''xb 

f H 
u uu 

f ^ 
u 

• 

V 

XX 

= «xx" 

f V 

X XX 

V. 

-<»UX  " C Huu'^(«ux  * C ''xx> 


-V  ^ f H 1'  ^ V 
Xt|.  u uu  u xt^ 


(64) 


where  all  quantities  are  evaluated  at  X»  and  n*  unless 

otherwise  specified,  The  terms  designated  with  the  bar  above 
indicate  the  values  along  the  nominal  trajectory  for  which 
that  variable  is  used.  The  optimized  control  for  the  next 
iteration  is  given  by 


u(t)  = u*(t)  + 3j^(t)  dx(t)  + 
+ tS^Ct)  dt^ 


(65) 


where  p^(t),  3^(t)  are  given  ast 

ei(t)  = (H^,^  * 


t32(t)  = -Hyy  fy  Vx^ 
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SjCt)  . 


(66) 


The  functions  with  the  subscripts  Indicate  the  partial  of 

that  function  with  respect  to  that  subscript.  Note  that  db 
* 

and  dt^  are  zero  except  during  the  restoration  phase,  and  the 
third,  fourth,  fifth,  sixth,  and  eighth  equations  of  (64)  are 
used  only  during  the  restoration  phase. 

The  Hamiltonian,  H,  is  defined  as 

H = L + f (6?) 

The  necessary  condition  for  an  optimum  control,  u*.  is  that 
a 0.  Taking  the  partial  derivative  of  the  Hamiltonian 
with  respect  to  u yields  the  following  u*i 

V * * RA 

where  is  defined  as  the  eleventh  component  of  the  costate 
equation  in  the  equation  (64),  etc.  Since  Ug  is  an  argu- 
ment of  the  sin  and  cos,  the  Ug*  equation  is  transcendental. 

A root  finding  subroutine  is  used  to  solve  for  from  the 
following  transcendental  equation! 

RB(u2*)  + cos(u2*)(Vq  1/2  p^e'^^T/^o^ 

• V,j,  a^m^  cos(Y,p))  + sin(u2*) 

^^ctT  ~ ® (^9) 

The  boundary  conditions  required  for  the  dif- 
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ferential  equations  (64)  aret 


a(t^) 


= 0 


F„  + T) 

•rV 


'xb 


xt. 


'bt. 


'bb 


XX 


Vf 


= '1^, 


«x  " ^xx  ^ 


f 


= 0 


= + B <r 


XX 


<Kx'  <f-  ''=« 


(70) 


where  < > signifies  the  inner  product,  and  b is  the  nominal 

Lagrange  multiplier.  The  computational  procedure  for  the 
program  is  outlined  in  the  next  chapter. 
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VI.  Computational  Procedure 


In  order  to  understand  more  clearly  the  computational 
procedure  used  for  the  optimal  missile  evasion  differential 
dynamic  programming  model  a synopsis  of  the  program  will  now 
be  discussed.  As  was  stated  earlier  in  this  thesis  the  main 
concern  of  this  program  is  to  determine  the  minimum  of  the 
missile,  as  well  as  the  optimum  controls  used  by  the  aircraft. 
The  Pj^  obtained  is  the  minimum  for  a particular  launch  condi- 
tion when  the  evading  target  uses  optimum  controls:  angle  of 
attack,  bank  angle,  and  coefficient  of  thrust.  In  the  real 
world  situation  the  Pj^,  for  the  same  launch  condition,  would 
be  equal  to  or  greater  than  that  Pj^  obtained  by  this  program 
depending  on  how  skillful  the  evading  target  pilot  was. 

For  this  program  the  terminal  condition  is  known,  that 
is,  the  stopping  criteria  occurs  when  the  LOS  from  the  mis- 
sile to  the  target  intercepts  the  missile's  fuzing  cone  as 
discussed  in  Chapter  IV,  The  duration  over  v/hich  the  control 
is  to  be  applied,  however,  is  not  fixed.  The  interval  [t^, 
t^]  is,  therefore,  not  specified  explicitly.  An  initial  time 
interval  of  one  second  was  selected  to  test  this  program.  A 
nominal  control,  u,  is  then  loaded  into  an  array  for  the 
three  controls,  Each  control  is  selected  as  a constant 
throughout  the  Interval  using  a step  size  of  .1  seconds, 
with  the  nominal  angle  of  attack,  u-j^,  of  zero  degrees,  the 
nominal  bank  angle,  u^,  of  zero  degrees,  and  the  nominal 
coefficient  of  thrust,  u^,  of  .025.  The  nominal  value  of 
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the  coefficient  of  thrust  is  based  upon  a 55,000  lb  airplane 
at  20,000  ft  using  military  power  (7000  lbs  of  thrust). 

A brief  outline  of  the  method  used  in  the  algorithm  is 
depicted  in  Pig.  8}  and  a more  detailed  analysis  of  the  pro- 
gram is  described  in  steps  1 through  6. 


Step  1 - Backward  Integration  of  the  State  Equations  Using 
Nominal  Control 

Since  it  is  known  approximately  where  the  terminal  con- 
ditions occur,  these  terminal  conditions  are  used  to  initial- 
ize the  state  equations  in  order  to  integrate  the  states  back- 
wards in  time  the  one  second  interval.  Referring  to  Fig.  9 
for  a geometrical  interpretation  of  the  terminal  conditions 
selected  for  this  problem,  the  initialization  of  the  state 
equations  for  the  backwards  integration  are* 


x„  - Xj  = 5 ft 

- Yt  = -8.66  ft 
Z„  = 20,000  ft 
= 20,000  ft 
Vj^  = 1866.4  ft/sec 
V,p  = 829.5  ft/sec 
ijfnj  = 30  deg 
=0  deg 

Ym  ^ ® 

Yt  =0 
(x^  =0  deg 

= .025 


(71) 
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Restoration 

Phase 


Bkwd  Integration 
(Constrained 
Terminal  Condition] 


Calcula- 
tion of 
db  & dt-f 


Step  Size/Linearization 
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^optimum 


Fwd  Integration 
of  States 
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Step  Size/Linearization 
Check 


Optimum  Cost 


Pig.  8.  Differential  Dynamic  Programming 
Model  Flow  Chart  (Concluded) 
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The  state  equations  are  then  integrated  backward  from  t^  to 
t^  using  the  nominal  controls  and  using  the  above  starting 
conditions.  The  backward  integration  is  stopped  at  t^  where 
values  for  the  states  are  obtained  and  designated  as  XI (t^) 

Step  2 - Forward  Integration  of  the  State  Equations 

The  state  equations  are  integrated  from  t^  to  t^  using 
the  same  nominal  controls  and  using  the  initial  conditions 
XI (tp)  as  determined  from  step  1.  During  this  forward  inte- 
gration the  nominal  cost  is  determined  using 

^non  ° V + ^ tf) 

♦ 1/2  A (RA-(Uj^)*+RB-(U2)^+RC.(Uj)^)dt  (72) 

t*' 

0 

^ere  F(X(t^)i  t^)  is  the  terminal  cost  as  discussed  in  Chap- 
ter III,  iif(X(t|.)i  tj)  is  the  terminal  constraint  as  discussed 
in  Chapter  IV  and  5 is  the  nominal  Lagrange  multiplier  which 
was  set  equal  to  1.  The  procedure  used  for  evaluating  the 
terminal  cost  and  constraint  is  discussed  in  step  3>  The 
states  are  then  evaluated  at  t^  and  designated  as  X(t^). 

Step  2 - Backward  Integration  with  Unconstrained  Terminal 
Condition 

The  following  set  of  differential  equations  are  inte- 
grated backwards  from  t^  to  t^i 

i 

States  I X 
A 
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(73) 


The  starting  conditions  for  integrating  the  above  equations 

A. 

backwards  are 


(74) 


The  terminal  conditions  for  the  costate  (V^)  equations  in- 
volve taking  the  first  partial  derivative  of  the  terminal 
cost  (Fjj)  and  the  terminal  constraint  (i|;^)f  while  the  termi- 
nal  conditions  for  the  equations  involve  taking  the  second 

partial  derivative  of  the  cost  {T  ) and  constraint  (i|f„„). 

Both  the  terminal  cost  and  terminal  constraint  are  functions 
of  a large  ntamber  of  the  states,  states  one  through  eleven 
for  the  terminal  cost,  and  states  one  through  ten  for  the 
terminal  constraint.  Since  both  the  first  and  second  partial 
derivatives  of  these  functions  are  required,  a large  number 
of  terms  are  involved.  Since  high  accuracy  is  not  required 
for  the  values  used  for  the  initialization  of  the  differen- 
tial Equations  (74),  the  analytic  solutions  to  the  initiali- 

• V 

zation  of  the  and  equations  will  not  be  attempted. 

The  derivative  is  approximated  using  divided  differences 


. F(Xo  4 AXp  - F(X^  - AXj. ) 
*1  2aX^ 


i = 1,2,  ...11 
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and 


+ aX^)  - F{Xq  - AX^) 
2aXj^ 


i = 1,  2,  ...10 


(75) 


where  F is  the  numerical  value  of  the  first  partial  deriva- 

X *tll 

tive  of  the  terminal  cost  with  respect  to  the  i state}  and 

ih  is  the  n\imerical  value  of  the  first  partial  derivative  of 

”Xi 

the  terminal  constraint  with  respect  to  the  i state.  The 
nominal  vector,  X^,  is  the  value  of  the  states  at  the  final 
time.  aX^  is  set  equal  to  1^  of  the  value  of  that  state  at 
the  final  time  plus  an  epsilon.  The  epsilon  is  added  in  or- 
der to  provide  a usable  perturbation  even  when  the  nominal 
value  is  small  or  zero.  This  prevents  division  by  zero.  The 
particular  values  for  the  arei 


AX^  s .01  X^(t^)  + 1 ft 

AXg  = .01  XgCtj)  + 1 ft 

AX^  = .01  X^Ct^)  + 1 f t 

A\  = .01  X^(t^)  + 1 ft 

AX^  - .01  X^(t^)  + 1 ft/sec 

AX^  = .01  X^(t^)  + 1 ft/sec 

AXy  a .01  X,p(tj)  + 1 degree 

AX0  = .01  Xg(tj,)  + 1 degree 

aX^  = .01  X^(tj)  + 1/2  degree 

AXio  = .01  Xj^Q(t^)  + 1/2  degree 

AXii  = .01  X^2^(t^)  + 1/2  degree  (76) 


Subroutine  FFX  calculates  F^  , and  subroutine  FPX  calculates 

^i 
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The  derivative  V is  approximated  using 

JvJv 

- P(X^  + 4X^  - 4Xj)  + F(X^  - iXi  - iXj)]  / 

4 AXi  AX . 

for  i / j i = 1,2,  ...11 
j = 1,2,  ...11 

and 

F(Xq  + AXj^)  - 2 P(Xq)  + F(Xo  - aX^) 


Vj 


AX, 


for  i = j i = 1,2,  . . ,11 


(77) 


where  F^  ^ is  the  numerical  value  of  the  second  partial  de- 

rivative  of  the  terminal  cost,  ilj^  „ , the  second  partial  de- 

i j 

rivative  of  the  terminal  constraint,  is  computed  in  the  same 

way.  The  aX.  and  AX.  are  the  same  values  as  that  previously 
^ J 

determined  in  Eq  (76).  Subroutine  FFXX  calculates  ^ and 
subroutine  PPXX  calculates  tlf^  ^ . 

Vj 

As  previously  stated  in  step  2,  at  the  end  of  the  forward 
integration  the  terminal  cost  and  terminal  constraint  values 
are  required.  For  the  terminal  cost,  subroutine  FFXX  is  used 
to  calculate  this  value.  The  terminal  cost  is  calculated  in 
that  portion  of  the  routine  where  zero  displacement  of  the 
states  is  required.  The  terminal  constraint  is  calculated  by 
subroutine  PPXX  using  the  same  procedures. 
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r 


I 


J j1|>  nil . 


Setuming  to  the  main  program  with  the  hoimdary  condi- 
tions, Equations  (73)  are  integrated  backward  from  t^  to  t^. 
Subroutine  F is  used  for  this  backward  integration.  During 
this  backward  integration  H is  minimized  with  respect  to  u 
to  obtain*\i*.  As  previously  stated,  the  equation  for  Ug  is 
transcendental,  therefore,  subroutine  F determines  U2*  from 
the  transcendental  equation  at  each  step.  The  solution  for 
U2*  is  obtained  by  iteration  using  the  root  finding  approxi- 
mation 


f(u2  ) + ^'(^2  ^ ^ ® 


or  solving  for  6U2 

f(u2  ) 


vdiere  6U2  is  the  change  in  U2  from  one  iteration  to  the  next 

and  f(uo  ) is  the  transcendental  fvinction  evaluated  at  u«  . 

o ^o 

f*(u9  ) is  the  derivative  with  respect  to  u,  evaluated  at 
^o  ^ 

U2  . The  new  updated  estimate  is  then 


U2  + 5U2  (80) 


This  iterative  type  procedure  continues  until  the  solution 
converges  (i.e.,  6U2  becomes  less  than  a specified  tolerance). 
If  |f*(u«  )1  approximately  equals  zero  then  the  second  deriv- 
ative  has  to  be  taken  and  6U2  becomes 


(81) 


6U 


2 


SS  + 


t/T^^ 

* f(u2  ) 


irtiere  the  proper  sign  is  chosen  such  that  ju^  + 6u^|  is  min- 
imized, 

• • 

The  and  equations  are  integrated  backwards  usir^ 
the  values  of  u*.  The  A equation  is  integrated  backward 
using  -ttie  difference  between  u*  and  Uj^oj^inal*  ^ 
eter  value  of  .0005  vfas  preselected  for  A.  Throughout  the 
entire  backward  integration  the  values  of  0Q^(t),  which  is 
defined  by  £q  (66) « are  stored.  The  time  is  noted  when  A 
exceeds  the  preset  value  of  .0005»  and  is  designated  as 
NEFF(2),  the  time  of  an  "effective"  change  in  A. 

Step  4 - Forward  Integration  of  State  Equations  Using  Improved 
Control 

The  state  equations  are  now  integrated  forward  again 
from  t^  to  t^  maintaining  the  same  initial  conditions  but 
using  the  new  control 

* u*  + dx(t)  (82) 

idiere  3j^(t)  was  calculated  in  step  3 and  dx(t)  is  the  differ- 
ence between  the  states  using  the  nominal  control  and  the 
states  using  the  control  designated  th^ring  the 
first  iteration,  dx(t^)  equals  zero  and  then  dx  is  calcu- 
lated for  the  remaining  steps  to  t^.  Throughout  the  entire 
forward  integration  the  new  coot,  designated  as  is 


calculated.  includes  the  integral  cost,  the  terminal 

cost,  and  the  terminal  constraint  penalty  calculated  the  same 
way  as  in  step  Z\  but  all  functions  are  evaluated  using  vi 
The  requirement  for  a i=>c.tisfactory  new  control  is  that  the 
change  in  cost  be  greater  that,  some  constant  times  A(t^). 
Jacobson  and  Mayne  Ref  (li25)  suggest  initially  that  the 
ccn.'«tant  equal  .51  therefore,  the  requirement  for  a valid 
new  control,  is 

Jnom  - Jstar  i <=  A<*o>  <83> 

If  inequality  (83)  is  satisfied  then  the  control  is 
loaded  into  the  H^jQ^inal  step  2 is  repeated, 

If  inequality  (83)  is  not  satisfied  then  a step  size  adjust- 
ment method  is  required  to  prevent  overstepping  the  region 
of  linearity,  The  step  size  adjustment  method  consists  of 
determining  the  position  on  the  integration  interval  half 
way  between  NEFP(2),  which  is  the  first  point  on  the  back- 
ward integration  whore  A exceeded  .OOO5,  and  t^.  This  new 
position  is  called  MK  in  the  main  program.  The  ^ control 
is  then  changed  by  the  follow -ng  method i 

^ew  * -nominal 
on  the  interval  from  t^  to  MK  and  then 

Hnew  “ (84) 

on  the  interval  from  MK  to  t|,.  Using  the  now  control,  step 
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4 Is  then  repeated  but  now  the  change  in  states,  dx(t),  Is 
aero  from  t^  to  MK  since  the  same  control  is  used  over  that 
interval,  and  on  the  interval  from  MK  to  tj,  dx(t)  is  cal- 
culated as  before.  Again  the  change  in  cost  is  evaluated 
and  the  criteria  for  a valid  new  control  is 

’^nom  - ’^star  ^ ^ ^(MK)  (85) 

The  step  size  adjustment  method  criteria  specifies  that 
the  change  in  cost  must  be  greater  than  a predetermined 
value.  If  the  criteria  is  not  met,  MK  must  be  moved  to- 
ward the  final  time  so  the  old  control  is  used  over  more 
of  the  total  interval.  This  iterative  type  procedure  is 
continued  until  MK  and  NEFP(2)  coincide.  When  this  condi- 
tion is  satisfied  then  C is  set  equal  to  zero.  When  C is 
set  equal  to  zei’o,  the  most  optimum  trajectory  has  been  found 
using  C equal  to  .5.  The  Uy^g^  control  is  then  again  loaded 
into  the  Bnominal  ^ repeated.  The  new 

criteria  for  a valid  control  is  that 


"^nom  ” *^star  ~ ® 


(86) 


This  is  a more  refined  criteria  in  that  as  long  as  is 

less  than  the  control  is  valid.  Again  this  iterative 
type  procedure  is  used  until  MK  and  NEFF (2)  coincide.  The 
most  optimum  -control  has  been  determined  by  treating  the  prob- 
lem as  a free  end  point,  that  is  the  terminal  constraint  cri- 


! /i. 


m 


(' 


c 


terla  has  not  been  satisfied.  The  next  step  will  be  to  re- 
store the  terminal  constraint  condition,  Eq  (55) • 

Step  ^ - Backward  Integration  with  Conetrained  Terminal 
Condition ^ 

The  following  set  of  differential  equations  are  inte- 
grated backwards  from  t^i  to  t^^i 

States  t X 

• 

V 

XX 

»xb 


'Vb 


f 


with  the  initial  conditions  ofi 


= Xn(tj) 

* F„  + '5 


XX 


'xb 


^XX  ^ ’^^xx 

= 11;/ 


^bb  = 


^xt^  “«x*Vxx^ 


(87) 
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+ 


(88) 


'Vt  - <J'x-  f>  * <f-  '^xx  f> 


The  initial  conditions  for  the  V^,  and  Vxb  equations  are 
determined  the  same  way  as  in  step  3.  The  initial  conditions 


for  the  and  ^ are  calculated  analytically  by 

subroutine  STCST.  During  the  backward  integration  using  Hnew 
32('t)  and  g^Ct)  are  determined  from  t^  to  t^  where  32('t)  and 
g^(t)  are  defined  byt 


S2<^)  = -«uu'^  C ''xb 

'"xt. 


(66) 


The  new  control  to  be  evaluated  on  the  next  forward  in- 
tegration is  given  by 


where  is  defined  in  Eq  (84)  and  where  db  and  dt^  are 

given  by 


(90) 


All  values  are  determined  at  t^,  therefore,  db  and  dt^  are 
constants  for  a particular  integration  cycle,  E is  also  a 
constant  for  a particular  integration  cycle  which  initially 
equals  1 as  suggested  in  Ref  (1j43).  The  quantity  dt^.  is  the 


change  in  the  final  time  of  the  trajectoryi  and  for  this 
problem  was  limited  to  a maximum  of  10^  of  the  final  time. 
That  is,  the  criteria  for  dt^  is 

Idt^l  < .1  tf  (91) 

where  t^  is  the  final  time  of  the  trajectory.  If  inequality 
(91)  is  not  satisfied  then  E is  set  to  E/2  and  then  db  and 
dtj  are  re-evaluated  by  Eq  (90).  If  inequality  (91)  is 
satisfied  then  proceed  on  to  step  6. 

Step  6 - Forward  Integration  of  State  Equations  Using  Optimum 
Control 

The  control  given  by  Eq  (89)  is  designated  as  u^p^  in 
the  program.  The  state  equations  are  then  integrated  from 
t^  to  tj  using  The  values  of  the  states  are  determined 

and  the  terminal  constraint  is  then  calculated  by  using  sub- 
routine FPXX.  The  value  of  this  optimum  constraint  is  des- 
ignated as  PCOSTP,  while  the  constraint  evaluated  in  step  4 
was  designated  as  PCOSTN,  In  order  to  determine  if  there  has 
been  an  improvement  in  the  end  point  error  the  following  cri- 
teria must  be  satisfied! 

1PC0STN1*>  1 PCOSTP 1 (92) 

where  again  PCOSTN  is  the  value  of  the  terminal  constraint 
by  using  Uj^g^»  and  PCOSTP  is  the  value  of  the  terminal  con- 
straint using  Uqp^.  If  inequality  (92)  is  not  satisfied  then 
set  E to  E/2  and  return  to  Eq  (90)  where  db  and  dt^  are  re- 
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evaluated.  If  inequality  (92)  is  satisfied  then  the 
control  is  loaded  into  the  Hj^ominal  3 

repeated.  The  stopping  criteria  for  this  restoration  phase 
is  when  the  terminal  constraint  satisfies 

|ilf(X(tj.),  t^)|  < .005  (93) 

The  final  values  obtained  by  this  program  are  the  fol- 
lowing I 

1.  The  trajectory  of  the  missile, 

2.  The  optimum  trajectory  of  the  evading  target, 

3.  The  optimum  controls  used  by  the  target  and, 

4.  The  value  of  the  terminal  cost,  which  is  the 

minimum  P,  the  missile  would  obtain  for  the 
k 

initial  conditions  as  outlined  by  Eq  (71). 

A detailed  schematic  of  the  differential  dynamic  algo- 
rithm used  for  this  program  is  illustrated  in  Fig.  10,  Ap- 
pendix A. 
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VII.  Results 


As  previously  stated  this  study  attempted  to  determine 
the  optimal  controls  to  he  used  by  an  aircraft  evading  an 
air-to-air  missile;  with  the  controls  consisting  of  the  angle 
of  attack,  the  bank  angle,  and  the  coefficient  of  thrust. 

The  state  vector  for  this  missile  evasion  problem  included 
the  distance  components  between  the  target  and  missile,  the 
velocity  of  both  target  and  missile,  the  heading  of  target 
and  missile,  the  flight  path  angle  of  target  and  missile,  the 
angle  of  attack  of  the  target,  and  the  coefficient  of  thrust 
of  the  target.  The  performance  of  the  system  was  measured  by 
the  minimization  of  the  terminal  cost  function,  which  was  the 
Pjj  of  th<=‘  missile. 

The  objective  of  the  terminal  cost  function  was  to  make 
the  model  correspond  to  detailed  simulation  of  the  end  game 
which  incorporated  the  shape,  orientation  and  vulnerability 
of  the  aircraft.  The  individual  vulnerability  of  a particular 
aircraft  was  not  considered,  but  only  the  generic  character. 
The  ellipsoid  parameters  were  chosen  to  fit  detailed  sim- 
ulation results  reasonably  well.  Previous  optimal  missile 
evasion  investigations  modeled  the  terminal  cost  function  as 
the  magnitude  of  miss  distance  only  and  did  not  consider  the 
vulnerability  of  the  aircraft  as  a function  of  the  orienta- 
tion geometry. 

Due  to  the  complexity  of  the  algorithm  for  this  optimi- 
zation program  and  the  imposed  time  limit,  the  desired  optimal 
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solutions  for  this  investigation  were  not  obtained. 

Prior  to  attempting  this  optimal  missile  evasion  problem 
the  following  system  consisting  of  two  states  was  solved  using 
both  the  first  order  and  second  order  methods, 

*2  “ "*1  - *1^  + « 

with  the  terminal  constraint  consisting  of 

tlf(X(t^)j  t^)  = X3^  - 2X3  - 2 (9^^) 

and  the  integral  cost  consisting  of 

L = 1/2  J + Xg^  + u^)dt  (95) 

Both  the  first  order  and  second  order  techniques  were 
attempted  in  order  to  evaluate  the  feasibility  of  using 
either  of  these  techniques  for  the  optimal  missile  evasion 
problem,  which  consists  of  a system  with  12  state  variables. 

In  this  sample  problem  it  was  determined  that  the  rate  of 
convergence  using  the  first  order  technique  was  extremely 
slow,  Witn  this  in  mind  it  was  determined  that  if  the  first 
order  method  was  attempted  on  the  evasion  problem  then  the 
slow  rate  of  convergence  would  make  it  extremely  difficult 
to  obtain  the  many  desired  solutions,  thus  the  second  order 
technique  was  selected  as  the  best  candidate  method  for  this 
program. 
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The  greatest  disadvantage  of  using  the  second  order  tech- 
nique is  the  high  dimension  of  the  system  of  differential 
equations.  For  any  optimization  problem  with  an  unspecified 
final  time,  using  the  technique  outlined  in  Ref  (li46),  the 
following  equations  must  he  calculated: 


-a 


H - H(X.  u,  t) 


-V. 


xh 


= (f^  * fu  Pi)''  V, 


'xb 


-V 


xt. 


-V. 


bt . 


-V  ^ f H f ^ V 

'^xb  ^u  "uu  ^u  ''xt. 


-V. 


bb 


-V, 


XX 


, _Y  T f ij  f ^ V 

''xb  ^u  ^uu  ^u  ''xb 

’ »xx  * fx’'  ’^xx  " '^xx  ^x 

"(*^UX  ^xx^  ^uu  ^^ux  ^xx^ 


-V 


tftf 


-V  f H f ^ V 

xtj  U UU  U Xtj. 


(64) 


with  the  control  determined  by» 

u(t)  = u*(t)  + 0j^(t)  dx(t)  + 

+ g^(t)  dt^  (65) 

and  where  P3  calculated  byi 

Pl^^)  = -«uu"^  <Hux  ^u""  \x) 

PgCt)  - fy 
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(66) 


C S3(t)  = 

In  Eq  (64)  the  total  number  of  differential  equations  to 
be  integrated  exceeded  I30  for  this  problem.  The  V„„  equations 
require  the  second  partial  derivatives  of  the  Hamiltonian  with 
respect  to  the  states,  and  several  of  these  equations  were  in 
excess  of  15  lines  of  fortran  computer  coding. 

In  addition  to  the  integration  of  the  above  equations, 
the  terminal  conditions  had  to  be  computed  for  each  equation. 
These  computations  required  the  use  of  a numerical  differ- 
entiation technique  as  was  outlined  in  Chapter  VI.  It  should 
be  noted  that  this  procedure  requires  a significant  amount  of 
calculation  in  order  to  obtain  the  numerical  derivatives  for 
• these  terminal  condition  values.  Along  with  these  calcula- 

tions, an  iteration  was  required  to  solve  the  transcendental 
equation  for  the  optimizing  control,  as  discussed  in  Chapter 
VI. 

When  attempting  to  run  this  prog.vjs  on  the  computer  the 
size  of  the  program  was  such  that  the  core  memory  allocation 
requirement  was  in  excess  of  177,000  (octal)  words.  The  cor- 
responding compilation  time  for  each  run  was  one  hundred 
thirty- two  seconds. 

In  order  to  reduce  the  execution  time  the  initial  flight 
inter-^ml  was  limited  to  one  second  with  the  intention  of  in- 
creasing the  limit  when  the  program  was  working  properly. 

The  results  to  this  point  indicate  that  the  state  equations 
L have  been  partially  successful  in  both  the  backward  and 
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forward  integration  for  the  one  second  time  interval.  How- 
ever, difficulty  was  encountered  during  the  backward  integra- 
tion of  the  costate  equations  which  determine  the  sensitivity 
coefficients,  The  backward  and  forward  integration  of  the 
state  equations  plus  the  initial  backward  integration  of  the 
costate  equations  required  over  fifty  seconds  execution  time. 

Due  to  the  large  memory  requirement  and  the  lengthy  exe- 
cution time,  this  program  was  placed  near  the  bottom  of  pri- 
orities in  the  central  computer.  Tum-around  time  for  each 
computer  run  varied  between  2 and  3 days. 

In  addition  to  the  slow  computer  tum-aroimd  time,  there 
exists  within  the  program  itself  the  inherent  problem  associ- 
ated with  debugging  the  complex  program  code.  That  is,  there 
exists  some  probability  of  programming  errors  even  though 
every  attempt  was  made  to  limit  such  a possibility.  The  com- 
plexity of  the  equations  used  in  this  program  definitely  con- 
tributes to  the  possibility  of  a mistake  in  derivation  of  the 
equations. 

This  investigation  demonstrates  the  desirability  of  find- 
ing an  alternate  method  of  genei'ating  the  equations  required 
for  the  second  order  optimization  technique.  The  alternate 
method  should  not  require  the  lengthy  derivation  of  partial 
derivative  equations  nor  the  large  amount  of  memory  required 
to  compute  them. 

A procedure  that  satisfies  these  requirements  is  a higher 
order  computer  language  designated  PROSE,  Ref  <2j1).  It  eval- 
uates automatically  partial  derivatives  as  a by-product  of 
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the  computations  in  the  model.  The  resulting  derivative  values 
are  exact  to  the  precision  of  the  computer,  as  though  they 
had  been  evaluated  by  formula  derivation  and  evaluation. 

This  derivative  evaluation  by  PROSE  would  eliminate  the 
lengthy,  time  consuming,  and  possibly  erroneous  derivation  of 
the  partial  derivatives. 
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I VIII.  Conclusions  and  Recommendations 
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The  method  of  solution  was  one  of  the  primary  concerns 
in  this  study.  The  rate  of  convergence  using  a first  order 
optimization  technique  for  solutions  to  this  missile  evasion 
problem  was  estimated  to  be  extremely  slow  and,  therefore, 
luiv^cceptable . The  second  order  technique  used  for  this  pro- 
gram also  has  disadvantages}  the  two  primary  ones  are  long 
computer  execution  times  and  the  large  ntimber  of  equations 
that  had  to  be  derived.  The  following  procedures  are  recom- 
mended as  additional  efforts  in  solution  methodology: 

1.  Simplify  the  problem  formulation  as  follows: 

a.  Use  an  altitude  difference  between  the 
missile  and  target  as  one  state  instead 
of  the  state  for  each  missile  altitude 
and  target  altitude, 

b.  Assume  the  time  responses  for  the  angle 
of  attack  of  the  target  and  the  coef- 
ficient of  thrust  of  the  target  are 
rapid  enough  to  neglect  any  time  delay 
in  the  response. 

With  the  above  simplifications  incorporated 
into  the  program  the  state  vector  would  be 
reduced  to  nine  states. 

2.  Complete  the  second  order  optimization  routine 
as  derived  in  this  thesis  in  order  to  obtain 
the  optimum  control  strategies. 
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J.  Investigate  use  of  the  PROSE  programming  lan- 
guage for  solving  the  evasion  problem  and  ob- 
taining the  optimum  control  strategies. 
Compare  the  results  obtained  with  the  optimi- 
zation program  derived  for  this  thesis  to 
those  obtained  with  PROSE. 
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Appendix  A 


Detailed  Schematic  of  the  Program 

A detailed  flow  chart  of  steps  1 through  6 used  in  the 
program  is  illustrated  in  Fig.  10.  The  initial  values  of 
the  variables  MK.  C,  b*  and  E are  annotated  at  the  top  of 
the  diagram.  Recall  that  the  initial  value  for  is 
given  as  the  following! 


\x^  ■ 0. 

Ug  * 0.  (96) 

u^  * .025 


where  the  controls  are  assumed  constant  throughout  the  en- 
tire time  of  flight. 

NEPF(2)  is  the  point  at  which  A exceeds  the  present 
value  of  .0005 I end  in  the  calculation  of  db  and  dt^  recall 
that  these  values  are  determined  by  the  following i 


-E 

-1 

V 

^t  t 

tftf^ 

s 

wm  m 

(90) 
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Pig.  10.  Detailed  Flow  Chart  of  the  Program 
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Pig,  10.  Detailed  Flow  Chart  of  the  Program 
(Continued) 
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Fig,  10,  Detailed  Flow  Chart  of  the  Program 
(Concluded) 
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